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Abstract 

It was first suggested by Englander et al to model the nonlinear dynamics of DNA 
PQ ' relevant to the transcription process in terms of a chain of coupled pendulums. In a 

related paper [4] we argued for the advantages of an extension of this approach based on 

considering a chain of double pendulums with certain characteristics. Here we study a 
1-^ . simplified model of this kind, focusing on its general features and nonlinear travelling 

j^r< wave excitations; in particular, we show that some of the degrees of freedom are 

actually slaved to others, allowing for an effective reduction of the relevant equations. 

(N' 

> 

(N 

Cp ■ In a seminal paper appeared a quarter century ago, Englander, Kallenbach, Heeger, 

Krumhansl and Litwin [12] suggested that DNA solitons, i.e. nonlinear mechanical ex- 

\Q \ citations of the DNA double helix, could have a key role in DNA functional processes - 

^£ ' such as DNA transcription and replication - in that they would provide automatic fo- 

O ' cusing of energy and synchronization of opening along the chain. They also were the 

r^ , first to suggest a simple mechanical model to illustrate their argument; this consisted of 

a one-dimensional chain of simple pendulums, each of them coupled to the neighboring 

ones. In this way travelling solitons were described by the sine-Gordon equation, which 

is known to support both topological and dynamical solitons. (Note that, contrary to the 

kS ' Davydov soliton in alpha-helices [10] which is intrinsically quantum, in this case we deal 

C^ ■ with classical mechanics and hence classical solitons.) 

Following the suggestion of Englander et al., a number of models for the nonlinear 
dynamics of DNA have been elaborated by different scientists in later years. 

Research was pursued essentially in two directions: on the one hand, in DNA denatu- 
ration the relevant aspect is the separation between the two helices, and one looks mainly 
at the degree of freedom describing distance between the two bases in a (Watson-Crick) 
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pair. Considerable progress in this direction has been obtained via the model introduced 
by Peyrard and Bishop [25]; this has then been refined by Dauxois [9] and extended in 
the BCP model by Barbi, Cocco and Peyrard [1, 2, 7, 23] (related models are discussed in 
[21, 35]). See e.g. [3, 24] for a discussion of results and recent advances in this direction, 
and [8, 26, 30] for matters related to (thermodynamic) stability and bubble formation. 

On the other hand, models for the roto-torsional dynamics of DNA have been considered 
by other authors; the torsional degrees of freedom are relevant in connection with the 
opening of the DNA double helix taking place to allow RNA-Polymerase to access the 
base sequence to transcript genetic information [34]. Several models have been proposed 
in this direction; references and some detail on these can be found in [19, 34].^ 

A particularly simple model (Y model) was proposed by Yakushevich [32, 33]; see also 
[13, 14, 16, 19]. Despite its simplicity, the Y-model succeeds in describing several relevant 
features of the DNA nonlinear (and linear) dynamics related to the relevant processes 
[19, 34, 36], and is thus the subject of continuing interest. In recent works it has been 
shown that all the (sometimes, very crude) approximations introduced by Yakushevich in 
her model - some of these substantially affecting the dispersion relations for the model 
[20] - have very little impact on the fully nonlinear dynamics and in particular on solitons' 
shape [17, 18]. 

The Y model should in any case be seen as only a first step in a hierarchy of increasingly 
accurate models [33], and has some drawbacks which should be removed by considering 
more detailed versions of the model; these are both quantitative and conceptual. 

On the quantitative side, we mention the impossibility to fit the observed speed of 
transversal waves in DNA with physical values of the parameters [36]; and the fact that 
soliton speed remains - provided it is smaller than a maximal speed - essentially a free 
parameter [15]. 

On the conceptual side, solitons are possible in these models thanks to the homogeneous 
character of the chain; but we know that actual DNA is strongly inhomogeneous, as 
bases are quite different from each other, and homogeneous DNA would not carry any 
information.^ 

The model we consider here describes the DNA double chain via two (rotational) de- 
grees of freedom per nucleotide, hence it will called a "composite Y model", following 
the nomenclature introduced in [4]. These degrees of freedom are related separately to 
rotation of the sugar-phosphate backbone (which can be of any magnitude) on the one 
hand, and of the nitrogen bases (which are constrained to a limited range) on the other 
hand. Both rotations are in the plane perpendicular to the double helix axis. 

This model is a simplified - and somehow a "skeleton" - version of the more realistic 
(and involved) model considered in [4]: in that paper we triggered our model to the actual 

^It should be stressed that these models deal with the DNA double helix alone, i.e. do not consider 
its interaction with the environment or other agents as RNA Polymerase; as such, they are not of direct 
relevance for processes involving other actors, albeit they are definitely a first step in the study of these 
(see e.g. [24, 34] for a discussion of the possible functional role of nonlinear excitations in DNA). On the 
other hand, the description they provide of the dynamics of the DNA molecule could nowadays be tested, 
in principles, by means of single- molecule experiments [22, 27]. 

^The idea behind considering such a model is of course that the homogeneous model can be considered 
as an "average" version of a more realistic model, which could be studied perturbatively; but in a model like 
the Y model (and more generally those considered in the literature) a perturbation breaking translational 
invariance would destroy the soliton solutions. 
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DNA geometry and dynamics, while here we consider a simplified model so to focus on the 
general abstract mechanism of interaction between topological and non topological degrees 
of freedom. For the same reason we will restrict to motions which are symmetric under 
the exchange of the two chains (as often done also in the analysis of the Peyrard-Bishop 
and the BCP models). 

The model we consider, and more generally the class of composite Y models [4], rep- 
resents an improvement with respect to the standard Y model both from a quantitative 
and qualitative point of view. 

On the quantitative, phenomenological, side, we have a much higher flexibility of the 
chain described by the model and dramatic consequences on the model ability to provide 
realistic physical quantities. E.g., with the composite Y model considered in [4], one 
obtains the experimentally observed transverse phonon speed using interaction energies of 
the physical order of magnitude^, and a selection of solitons' speed. 

On the qualitative side, which is maybe even more interesting (and more widely ap- 
plicable than merely DNA), the composite Y model is remarkable in that the uniform 
and the non-uniform parts of the DNA molecule (backbone and bases respectively) are 
described by separate degrees of freedom. It happens that, as a consequence of the ge- 
ometry of the DNA molecule, the degree of freedom describing the backbone supports 
topological - hence strongly stable - solitons, while the one describing the motion of bases 
performs quite limited excursion due to steric hindrances; as mentioned above, this is taken 
into account in our model. It is thus quite conceivable that introducing in the model a 
non-uniformity which affects only the latter degree of freedom, a perturbation approach 
would allow to obtain solutions in terms of perturbed solutions for the uniform model^; 
see sections 4 and 5 below. 

The perturbative approach is also attractive in that by a suitable limiting procedure, 
see section 3 below (and [4]), the composite Y model reduces to a standard Y model, 
for which exact solutions can be obtained. Thus, a perturbative description for the uni- 
form composite Y model can be obtained by perturbing the system near the standard Y 
solutions. 

A drawback of general composite Y models is that the equations describing its dynamics 
are too complex to be solved analytically, and even the perturbative expansion around the 
solitonic Y solution can be very hard to control [4]. On the other hand, most of the nice 
features of the composite model seem to be quite generic, related mostly to doubling of 
the degrees of freedom and their different topological features, i.e. largely independent 
of the model details. Numerical investigations for the "realistic" composite Y model of 
[4] showed that introducing a number of simplifications into the model does not affect its 
main features - confirming the observations for the standard Y model [17, 18] - but surely 
makes it easier to handle it at the analytical level. 

Motivated by the previous arguments, we study in this paper a simple - possibly the 
simpler - composite Y model, for which the comparison with the standard Y model is 
immediate. Our main purpose is indeed to focus on the essential features introduced by 

''As mentioned above, within the standard Y model the same speed can be fitted only using an intrapair 
interaction energy which is about 6000 times the physical value [36]. 

""It should be mentioned that the BCP model [1, 2, 24] also presents the interaction of topological and 
non-topological degrees of freedom; however, there the topological degree of freedom is actually a cyclic 
variable and one has correspondingly a conservation law: this leads to a dynamics less rich topologically. 
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Figure 1. Notation for cacli of the double pendulums along the two chains; see text. 

having two different kinds of degrees of freedom in chain models. 

For the sake of concreteness we make reference to - and discuss the consequences for - 
DNA dynamics, but it will be quite clear that most of our discussion is more general. 

1 The model 

Let us now describe our model. We define this in abstract terms, but it can be useful - in 
order to fix ideas ~ to recall that in applying it to DNA the "first pendulums" mentioned 
below and associated to angles 9 represent the nucleosides (segment of the phosphodiester 
chain and the attached sugar ring), while the "second pendulums" associated to angles (p 
represent the nitrogen bases. 

We consider two chains of double pendulums. On each chain axis there are equally 
spaced sites at points of coordinate z = nS, with 6 a dimensional constant and n G Z. 

At each site i G Z on each chain axis, there are identical simple^ "first pendulums" of 
length R and mass M, which can rotate with no limitations in the plane orthogonal to 
the axis. The rotation angle of the first pendulum on chain a = ±1 at site i will be 9'^. 

Attached to the point mass of each "first pendulum" there is a simple "second pendu- 
lum" of length r and mass m, which can rotate in the same plane as the first pendulums. 
The angle of rotation of the second pendulum (with respect to the direction of the first 
pendulum) will be (j)f. Details are shown pictorially in fig.l. 

The second pendulum is not free to swingle through a full circle, but is instead con- 
strained to stay in the range \(j)f\ < (pQ. This constraint will also be modelled by adding a 
constraining potential, i.e. a potential 14(0) which has the effect of limiting de facto the 
excursion of the (j) angles. 

The reason for this limitation in the </> range is our intention to model the DNA molecule: 
in that the bases cannot move outside a certain range, as they would otherwise collide 
with other parts of the molecule (the phosphodiester chain or the sugar ring) . 

®That is, pendulums consisting of a perfectly rigid and massless bar rotating about an extremum, and 
of a point mass fixed on the opposite extremum. 
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As for the relative position of the two double pendulums chains, this is such that the 
distance between the suspension points at the same site - and thus with the same z 
coordinate - is lA. Equivalently, at rest (i.e. for Q^ = (p'"' = for all a = 1,2 and 
i G Z), the distance between masses at the end of second pendulums at the same site on 
the two chains is Iq := 2(A — R — r). 

Let us denote, dropping for a moment the index i and writing a = it instead of a = 1, 2, 
as (x, y) the cartesian coordinates - in the plane orthogonal to the double helix axis z, 
with origin in the central point - of the the point mass at the end of the second pendulum. 

In the equilibrium position, given by = (?!> = 0, we have (x, y) = (xo,yo) with Xq = 
^A ± R±r, yo = 0. In general, 

x^ = =fA±Rcos9±rcos{9 + (j)) , y^ = ±Rsm9 ±rsm{9 + cp) . (1.1) 

Remark 1. In the standard Yakushevich model, one considers the approximation Iq -^ 0; 
it is known that this approximation has a strong impact on linearized dynamics and dis- 
persion relations [20], but for what concerns fully nonlinear dynamics it has been observed 
that - at least in the framework of the Yakushevich model - the approximation has a very 
little impact on travelling solitons [17]. This suggests the possibility of considering the 
same approximation Eq = 0, hence A = R + r, here as well. We will indeed adopt this. 

1.1 The Lagrangian 

Let us now describe the different terms appearing in the Lagrangian C for the two chains 
of double pendulums. 

The kinetic energy of the double pendulum (i.e. a disk and the attached pendulum) at 
site i on chain a is Tf = (1/2) [/(6'f)^ + m{{xf)^ + (yf)^)]; using (1.1) to express x" and 
yf in terms of 9 = 9'j^ and (/) = (pf, we get with simple algebra 

i;" = {1/2)9^ + {m/2)\R^9^ + 2Rr cos (/)( {9 f + 9<i))+r^ {9 + <i)f] . (1.2) 



The total kinetic energy is of course T = ^^ ■ T". 

We now pass to consider coupling between different double pendulums; we will have 
only nearest-neighbor interactions, which will be of three different types^. 

(1) A torsional coupling between successive disks on each chain, described by a potential 
Ut'^ = Kt[l-cosi9t+,-9'^)]; (1.3) 

the total torsional potential is of course Vt = '^Zai^t'^- This, albeit entirely natural in 
terms of the mechanical model (and thus mentioned here), is actually of lesser relevance 
in modelling DNA and will thus be overlooked; see Remark 3 below. 

'^In analyzing small amplitude dynamics, the "helicoidal" coupling between disks h sites away from 
each other (with 2h5 the pitch of the DNA helix) would also play a role [9, 13, 16, 19). Physically this 
interaction is mediated by Bernal- Fowler filaments [10], and it can be described by a potential 

[/« = (1/2) Kn [{e[li - e[-^f + {e\ii - ef'f] ; 

however, this is unessential in the fully nonlinear, i.e. large amplitude, regime and will hence not be 
considered here, as we focus on fully nonlinear excitations (in the small amplitude regime, introducing this 
term removes a degeneration of the system and is therefore qualitatively relevant) . 
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(2) A stacking coupling between successive pendulums on each chain, described by a 
potential Us = {l/2)Ks[{xf^i — xf)'^ + {yf^i — yf)"^]- This is expressed in terms of the 
9f and (j)f angles, using again (1.1) and with A = r/R, as 

^i") . iK.RV2) [(CO..? - c„,«f„ . A|co.„f + .?) - cos(,?, + 0^.^ 

+ (sin^f - sin^f+i + A[sin(0^ + Of) - sin(,^f_^i + ^f+i)]) ; 

the total stacking potential is of course Vs = ^^ j Vs'^ . 

(3) A pairing coupling between double pendulums at the same site on opposite chains. 
We assume this is due to an interaction between the point particles at the extremal points 
of the pendulums, described by a potential 

t/W = Fip,) (1.5) 

where pi is the distance between point masses of the two pendulums at site i on the two 
chains, /9j := [(x^ — x\ )'^ + {y\ Vi )^]^'^. The total pairing potential will of course 
heV, = Y.^ui^^. 

Recent work of ours in the framework of the simple Yakushevich model [17, 18] has 
shown that the shape of travelling topological solitons is very little affected by the detailed 
shape of the intrapair potential [18]. We will thus adopt the simplest choice, i.e. an 
harmonic potential Up = {l/2)Kppf; this becomes anharmonic when expressed in terms 
of the Of, (j)f rotation angles, 

[/« =AKp[{R + rf -rR{l- cos (l)'i)-{r + R) {rcos{(t)'i + Of) +Rcos Of)] . (1.6) 

Putting all these together, and considering also the constraining potential Vc to be 
discussed below, the Lagrangian describing our system will therefore be 

£ = T - (Vt + Vs + Vp + Vc) . (1.7) 

Some remarks are in order concerning the potential interactions introduced above. 

Remark 2. We have explicitly assumed harmonic potentials for the torsional and stacking 
coupling between successive elements (disks and pendulums respectively) on the same 
chain. It should be mentioned, in this context, that recent work by Saccomandi and Sgura 
has shown that the choice of a harmonic potential for the interpair interactions (torsional 
and stacking couplings) is not without effect: considering non-harmonic couplings leads to 
a stronger energy localization [28] . We will defer consideration of the effect of anharmonic 
interpair couplings in our model to later work. 

Remark 3. We also note that we have considered here both a torsional coupling origi- 
nating in the coupling between adjacent first pendulums (in DNA, in the phosphodiester 
chain) and a rotational interaction originating in interaction between adjacent second 
pendulums (in DNA this models the stacking interaction between neighboring bases). 

The physical origin of these in actual DNA molecule is rather diverse. Indeed, the 
stacking interaction is essentially due to polar interaction between the bases; these are 
essentially flat and rigid assemblies of atoms with a strong dipolar momentum, and cor- 
respondingly we have (in physico-chemical notation) a it — it bond. On the other hand. 
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the segment of a hypothetical phosphodiester chain without attached bases would be es- 
sentially free to rotate with respect to each other (if not for rather small polar forces), the 
main obstacle to such motion originating in hindrances caused by bases and their occupied 
volume and polar moments. 

It is thus reasonable, as announced above (and also in view of the considerable com- 
plications of equations that would be obtained from the full Lagrangian and the fact the 
two motions are strongly related geometrically), to consider the approximation Kf ^ 0, 
as we do from now on. This corresponds to the well known fact that the forces stabilizing 
the DNA double chain are essentially the interbase pairing interaction and the stacking 
interaction [24], i.e. those corresponding to our terms Vp and Vg- 

Remark 4. Further support to this kind of approximation comes from the results of [4], 
where the model with non vanishing torsional interaction has been investigated in detail. 
In that case, the torsional energy is a full order of magnitude lesser then the stacking 
energy. Moreover, numerical soliton solutions studied in that context showed that the 
form of the soliton remains almost the same if Kt = and all the torsional potential 
energy is transferred on the stacking term (see the Fig. 8 of [4]).* 

Remark 5. We stress that, as quite common in DNA mechanical models, we are dealing 
with the DNA molecule per se, i.e. we do not consider interactions of the double chain 
with the solvent (the fluid environment it is immersed in). This interaction would lead to 
exchange of energy with the solvent, with dissipation and random terms appearing in our 
dynamical equations (some authors have taken into account interaction with the solvent 
by introducing correction terms in the intrapair potential Vp [37, 38]). We hope a model 
modified in this direction can be studied in the near future. 

1.2 Constraints and constraining potential 

It should be recalled that the angles 6 and (j) are not on the same footing. Indeed, while 
9 is free to take any value, 6* G R, the angle 4> is constrained to be in the range 

<P G [-<Po,(l>o] ■ (1.8) 

The Euler-Lagrange equations (see below) for the Lagrangian (1.7) should be comple- 
mented with such a constraint. 

Note that (1.8) is an anholonomic constraint, and taking it into account is somehow 
inconvenient. It may be convenient, in particular in numerical computations, to remove the 
constraint and introduce instead a constraining potential. That is, we allow in principles 
(j) to take any value, but introduce a potential Uc{(f>) which makes very costly energetically 
any value of 4> outside the range specified in (1.8), and has an irrelevant effect on the 
dynamics when (j) is within this range. 

Any expression for the potential Uc would be acceptable, provided Uc is essentially flat 
for 1 01 < 1^0 and rises sharply as soon as we approach the border - or are outside - of this 
range. A convenient form for Uc is 

Uc{4>) = eM{Kc/2){<p' - cPl)] (1.9) 

*At the kinematical level the only difference between the model considered here and that in [4] is that 
here the base is a point particle whereas in [4] it is described by a disk of non vanishing radius. 
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with Kc >> 1. The total constraining potential appearing in the Lagrangian will then of 
course be 



1.3 Equations of motion 

The equations of motion for the model are the Euler-Lagrange equations for the Lagrangian 
(1.7). These are obtained with standard algebra, and are relatively involved expressions, 
which we do not write explicitly. The equation for the dynamics of variables associated 
with site i involve also values of variables at neighboring sites i ±1. 

We are mainly interested in solutions which vary slowly in space on the space scale set by 
the intersite distance 6. It is thus convenient to pass to the continuum approximation; this 
means promoting the array of values </>f (t) to a field variable ^"(z, t) such that ^""{nS, t) = 
(pnit), and similarly promoting 9^{t) to a field variable 0"(z, t) such that @"-{n5, t) = 0'^{t). 
(Note we are keeping the dimensional constant 5, so to have physical units - rather chain 
units - for z.) 

Moreover, in order to focus on the essential features of the model - and similarly to what 
is traditionally done in studying the Peyrard-Bishop and Barbi-Cocco-Peyrard models 
[1, 2, 24, 25] ~ we will just consider symmetric motions; that is, we assume (f>\ (t) = 
(p- (t), dl (t) = 0^ (t) which of course entails equality of the field variables as well. 
We will thus write simply ^^^>{z,t) = ^{z,t), Q^^'{z,t) = Q{z,t); and similarly for 
derivatives. 

Remark 6. It should be stressed that if we restrict to symmetric solutions in discussing 
the standard Y model, then only one type of kink solitons can exist, i.e. those connecting 
equilibrium positions at x = itcxD with a total phase shift of ib27r; higher solitons, even 
if with symmetric topological numbers, would require a deviation from a fully symmetric 
dynamics. 

In this setting, we write (j)i ^ ^ and 9i ~ Q for variables at site i; and for neighboring 
sites we have, expanding in Taylor series up to order two in 5, 

<Pi±i = ^±6^, + {i/2)d^^,, , ei±i = Q±6e, + {i/2)5^e,,. (i.ii) 

With these, the Euler-Lagrange equations issued from (1.7) are rewritten as second 
order PDEs for ^ and Q. In order to simplify the writing, we introduce the notations 

a(<I>) := r^ + rii cos ($) ; /?($) := r'^ + R^ + 2rRcos{^) . (1.12) 

The Euler-Lagrange equations are then given explicitly by 

d'^Ksir^ + rR cos <P) <^ccx + S^Ksir"^ + R^ + 2rRcos<^) @^^ + 

- [MR'^ +m{r'^ + R'^ + 2rR cos ^)]Qtt - m{r'^ + rRcos^)^tt = 
= ri?[(,52i^, sin $)$,($, + 26,) - (m$i($i + 29*)] + 

+ AKp{rR + R^) sin 6 + AKp{rR + r^) sin($ + G) ; , . 

S'^KsV^^^^ + 5^Ks{r^ + rRcos^)Q^^+ ^^'^''^ 

— TTir^ <5tt — m(r^ + ri? cos <I>) ©it = 

= - {rRsm^){AKp - rnQj + S'^K^Ql) + 

+ AKpr{r + R) sm{^ + Q) + iC, exp[(i^2/2)($2 - <^g)] $ 
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2 Travelling wave solutions (solitons) 

We are specially interested in travelling wave solutions, i.e. solutions satisfying 

cD(z,t) = y,{z-vt) = ifiO , &{z,t) = i){z-vt) = ^(0 . (2.1) 

These will be called "solitons" - provided they satisfy suitable boundary conditions, see 
below - following the common language in DNA modelling, and are indexed by an integer 
(i.e. a topological winding number). We stress that they are not proven to be dynamical 
solitons in the rigorous mathematical sense [6]; on the other hand, they are topological 
solitons in the sense of field theory [11]. 

2.1 Equations of motion 

With this ansatz, the equations (1.13) read 

ix{r'^ + rRcos(t)) (f)" + (J + 2^ri?cos0) 6*" = 

li{rRsm(l))[{(l)'f + 2(l)'e'\ - 4Kp[{rR + R^) sin 9 + {rR + r^)sm{(j) + 0)] ; 
fir'^ (j)" + fj,{r'^ + rRcos4>)0" = {4:KprR + firR{0')'^) smcp + 

- 4Kpr{r + R) sm{<f> + 0) + K, exp[{Kj2){<P^ -<Pl)]4>- 



(2.2) 



Here we have simplified the writing by introducing the parameters 

^ ■= mv^ - Ks6^ , J = {fi{r^ + R^) + MR^v^) . (2.3) 

In the analysis of the (2.2), we will proceed perturbatively; see next section. 

2.2 Finite energy condition and boundary conditions 

The equations (2.2) are a reduction of (1.13), which are in turn issued by the Lagrangian 
(1.7). The latter is well defined on fields ^{z,t),Q{z,t) satisfying a finite energy condition; 
in the continuum approximation the action S is given by integration over z £ (— oo, +oo) of 
the Lagrangian density obtained as the continuum limit of (1.7). In the present case, where 
we have a kinetic and a potential energy parts, the finite energy condition corresponds to 
requiring that for large \z\ the kinetic energy vanishes and the configuration correspond 
to points of minimum for the potential energy. 

By the explicit expression of the kinetic energy and of our potentials, this means 

<^{±oo,t) = , e(±oo,t) = 2n±7r , 

$,(±oo,t)=0, e,(±oo,t)=0, (2.4) 

$t(±oo,t) = 0, ei(±cx),t) = 0; 

where of course ^{±00, t) stands for lim2_»-|-oo ^iz,t), and so on. 

When we pass to the travelling waves equation, this means that we must impose on 
the functions (p{S,), 'd{S,) the boundary conditions 

ip{±oo) = , i?(±cx)) = 2n±7r , (p'{±oo) = , '(?'(±c5o) = . (2.5) 

In the following we will describe the dynamics of i^iO^V^iO i^ terms of motions (in 
the fictitious time ^) in an effective potential; note that such a motion can satisfy the 
boundary conditions (2.5) only if (99, t?) = (0, 27r/c) is a point of maximum for the effective 
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potential. The solutions satisfying (2.5) can hence be classified by the winding number 
n := n+ — n_. 

In the following we will in particular attempt to describe the basic (n = 1) soliton; note 
that - as stressed in Remark 6 - this is the only case for which we can fully compare the 
solutions of our model with those of the standard Y model within the frame of symmetric 
solutions. 

3 Perturbative expansion around Yakushevich solitons 

As remarked above, while (9 and hence Q and) -d can take any value, the range of (|0| and 
hence |$| and) \(p\ is limited by (po; thus if 0o = e << 1, we are guaranteed if will also be 
of order e. 

On the other hand, for (pQ ^ and hence (/? ^ 0, the (p degree of freedom is frozen 
and we are effectively back to consider the Y model [4]; we would then obtain the soliton 
solutions to our model as a perturbation of the corresponding soliton solutions for the Y 
model. The Y model can be also recovered from our composite model by acting on the 
geometrical parameter characterizing it [4]; this essentially amounts to a suitable limit 
r = 0. Here "suitable" means that some care should be taken in the rescaling in order to 
avoid a singular limit arises. 

We want to study the equations of motion for our model in the regime of small e, 
by looking at them as a perturbation of the standard Yakushevich equations. We stress 
that we are considering the symmetric (in the chain exchange) setting, and adopting the 
contact approximation A = R + r. 

We will take (po = e « 1 and expand in e. We will also take r = erg and, in order to 
keep the total length of the double pendulum A = R + r constant, R = A — erg. 

Remark 7. As shown in [17], one can obtain analytic results also without letting Iq -^ 0; 
on the other hand it is also shown there that a nonzero io does actually produce a very 
small modification in the solitons' shape, but yields much more involved analytic formulas. 
Thus, in the present context it is convenient to adopt Iq = 0. Note that in any case we 
should keep the physical distance A fixed as e is varied. 

We will expand 99 and 1? to order e. We write^ 

ip = 'ip + eipo +0{e^) , 



{} = 'dQ + erj + 0{e 



2\ 



(3.1) 



We plug this into equation (2.2) and expand the equations so obtained to first order in 
e, obtaining the equations 

A^{{pi - Mv^)i}o" - AKpsiwdo) + A [-(^ - Mv^){2roi}o" - Ar]")+ 
+ro(V'" + 2/U'!9o") cos V' — ip'ro{tp' + 2'i}Q')fj,simp — AKpr^ sim?o+ 
-4irp(ylr/cosi?o-2rosini?o)-4-ftrp^osin(V' + t?o)]e = 0; (3.2) 

-exp[(K,/2)V'2]i^e^ + [A^ro(cosV')V + ^ro(4irp + ^(t?o')^)smV+ 
-AAKpr^ sin(^ + i?o) - eM{Kc/2)il)^]K^{l + KcV'^)¥'o] e = . 
®The reader might me puzzled, as the title of this section seems in contradiction with the expansion 
(3.1): Yakushevich soliton involves only the topological angle i9, and here we are not even assuming tJo 
is the Yakushevich soliton. We will soon see that i/) = and i?o is precisely the Yakushevich soliton as a 
consequence of the equations of motion, i.e. (3.1) leads to an expansion around the Yakushevich soliton 
with no need for a specific assumption. 
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At order zero the second equation in (3.2) yields simply 

ip = . (3.3) 

With this, the equations (3.2) become quite simpler, i.e. reduce to 

A^dfi- Mv^)^o" - 4i^psim?o) + 

+A [2Mrov'^'do" + A{fi - Mv'^)r]" - AAr]Kp cos 'do] e = ; (3.4) 

[-Kcifo + Anroi^o" -AAKpro sin ^o] e = 0. 

We will now compute the leading term contributions for i? and for f, this will suffice 
to show our main point, emphasized in section 4. See Appendix B for the 0(e) correction 
to i}. 

3.1 The leading term for i!) 

As for the first of (3.4), at order zero we get a sine-Gordon equation, 

^o" = -{4:KpR^/J) sirn^o • (3.5) 

This can be seen as describing the motion (in the fictitious time ^) of a point particle of 
unit mass in the field of an effective potential 

Vo = ^— (l-cosT?o), (3.6) 

where the additive constant has been chosen so that for tJq = we get Vq = 0. 

It is important to stress that the leading term of the perturbative expansion reproduces 
the Y soliton equation (3.5) without any further constraint. In particular there is no 
constraint on the travelling speed of the soliton. 

Remark 8. This behavior should be compared with that described in [4] for a similar 
model, which differs from that under consideration here, essentially for the presence of a 
non vanishing torsional potential. In this latter case setting cp = allows to recover the 
Y soliton but the speed of the travelling wave is fixed in terms of the torsional constant 
Kf] moreover, the solitonic solution exists only when Kt/I < Kg/ra [4]. The fact that 
the above constraints disappear when Kt = Q sheds light on the physical meaning of the 
result of [4]: the constraint fixing the soliton speed found in [4] has a purely dynamical 
origin and is completely independent from both the geometry and the doubling of the 
degrees of freedom introduced in the composite model; its origin has to be traced back 
to the interplay between torsional and stacking potential energy which characterizes the 
composite Y model analyzed in [4]. 

As discussed above, the boundary conditions (2.5) require that i?o = corresponds to 
a maximum for the effective potential; it is apparent from (3.6) that this is the case if and 
only if J < 0, which we assume from now on. By the expression of J, see (2.3), this is the 
case provided 

"' < ^^^' (mI tm^ 2 (3-7) 
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Figure 2. Left (a): the function ^oiO^ see (3.9). Riglit (b): tlie function (fioiO^ see (3.11). Here 
we adopted the parameter values discussed in Appendix A, and plotted the case w = 0; we also 
used Kc = 30 in (b). The distance along DNA axis £, is measured in units of the intersite distance 
6, and do, ipo in units of tt. 



(the limiting speed for the standard Y model is recovered for r ^ 0). It should be 
noted that for v satisfying (3.7), the parameter n is also negative, /U < (more precisely, 
-Ks6^ <fi< -Ks5^{l + (m/M)(l + {r/Rf))-^.) 

The "conservation of energy" in the potential Vq for motions satisfying (2.5) reads 
(1/2)(t?^)2 + yo(i?o) = 0, i.e. we get t?o' = AR^{Kp/\J\) sin(T9o/2). The equation is 
integrated explicitly by separation of variables, yielding 



[log(tan(^o/4))] = i^(e - Co) 



(3.8) 



where we have defined K := 2R^JKp/\J\] note that choosing ^^o = corresponds to 
requiring that '!9o(0) = vr. With this, and inverting (3.8), we have the solution (plotted in 
fig.2.a) 



t?Q(^) = 4 arctan [exp (i^,f )] . 



(3.9) 



3.2 The leading term for Lp 

Let us now pass to consider the 0{e) terms in (3.4). In particular, the second equation 
reads 

Kcipo + 2qR{AKp sim?o + ^^o") (3.10) 

and yields directly an expression for ipi in terms of 7?o- Using the equation (3.5), this is 

SKpqRijiR'^ - J) 



^0 



JKr 



sim?o • 



(3.11) 



This function is plotted in fig.2.b. 
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4 Slaved fields and realistic DNA modelling 

We have thus obtained the first order term of the e expansion for 99 in terms of the t?o 
solution, with an algebraic rather than differential equation, see (3.11). In other words, 
ifQ is slaved to t?o; the same could be said more generally for the field <1> with respect to 
the field G. 

It is quite apparent that this feature is due to the different nature of the fields in 
our model. As already pointed out, is free to change and go round the circle (is 
topological), while $ is constrained to small variations. When we look for finite energy 
solutions, moreover, Q will join different vacua at z = ±00, while <I> can at most have 
small oscillations around the zero value. 

If we look at this feature from the point of view of DNA modelling, one point should 
be stressed. That is, suppose we had a non-homogeneous model, the inhomogeneity being 
however such to affect only the dynamical parameters related to the interaction between 
(j)i angles. In this case, the constants appearing in (3.11) would not be constant any more, 
and actually depend on the z coordinate - in the travelling wave reduction, on ^ and also 
on the parameter v describing the speed of the travelling wave. This would lead to a more 
involved expression for ipQ in terms of t9o) now depending explicitly on ^ as well as on v, 
and with all parameters being a function of ^; but the relevant point is that we would still 
be able to express 930 in terms of t?o- 

Needless to say, this is only true at first order, and higher order computations would 
soon become very much involved. On the other hand, it remains true that the "master" 
dynamics would be the one related to the 9 angles, and the (p angles would remain slaved 
at all orders; that is, writing 



00 



T? = i?o + e X] ^^^k ' '^ = X] ^^^^ ' 



fe=0 A:=0 

the equation for rf' would depend on ipo, ...,^fc_i (beside depending of course on t9o and 
on rji, ...,7]k-i), while ipk would be determined algebraically by i'&Ojrji, ...,7]k-i) in the 99 
expansion. 

The separation of the degrees of freedom in "master" and "slaves" is related to a 
nonperturbative feature of the dynamical system (2.2). The first integral of Eqs. 2.2) can 
be written in the form cp' = F{9' , 6, 93, E) where E is an integration constant ( the energy). 
Eliminating ip" from the system (2.2) and using the first integral above, one can always 
write the resulting equation in the form G{9" ,6' ,9,(p,E) = 0, which defines in implicit 
form 99 as a function of {9",9',9,E). Thus, once a (perturbative or non perturbative) 
solution for 9 is known, the corresponding expression for ip can be obtained algebraically. 

5 Discussion 

We have considered a mechanical model consisting of two chains of double pendulums, 
coupled with each other at each chain site. This "skeleton" model is the simpler specimen 
of a class of models for DNA roto-torsion dynamics which is aimed at improving the Yaku- 
shevich model, within the approach set forth by Englander et al. [12], and should show the 
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essential characteristics for the class keeping to a minimum the unessential mathematical 
complications. 

An element of this class - based on a geometrically more complex and realistic descrip- 
tion of the internal structure of DNA - has been considered in [4]. We have introduced 
two simplifications with respect to that model, which have been suggested by numerical 
results for that model itself: bases have been considered as point particles, and torsional 
interaction disregarded. These two quite reasonable approximations simplify drastically 
the dynamics of the system. 

We have determined the PDEs describing the motions in the long wavelength approx- 
imation, and also the reduction to travelling wave solutions. The finite energy condition 
imposes boundary conditions such that the finite energy solutions can be classified by a 
topological integer - actually a winding number for the field and hence for t? - and 
minimal energy solutions in each topological sector are topological solitons. 

The equation of motion describing travelling waves for our skeleton model are simple 
enough so that a (first order) perturbative expansion around the well known Yakushevich 
soliton can be used to find approximate analytical solutions. 

It seems appropriate to point out some drawbacks of our discussion. 

(A) The question if our model equations support dynamical solitons in rigorous mathe- 
matical sense [6] was not solved nor tackled. 

(B) Similarly, we did not investigate the convergence of the perturbation expansion around 
Y solitons, and just considered first orders contributions. 

(C) Also, as mentioned in Remark 5 and relevant for DNA modelling, we did not consider 
interactions of the double chain with the solvent. 

On the other hand, we provided a simple model for a nonlinear (double) chain showing 
some features which - in our opinion - are quite interesting both per se (i.e. in the frame 
of nonlinear dynamics for discrete systems) and for applications, in particular to DNA 
roto-torsional dynamics. 

In particular, the model is simple enough to allow for an explicit analysis, yet displays 
a very interesting interaction of topological and non-topological degrees of freedom. 

We have shown that the equations can be tackled in terms of a perturbation expansion 
around the Yakushevich soliton, and computed the behavior of the new (non-topological) 
field at first order. Quite remarkably, this is completely determined by the behavior of the 
topological field at order zero. That is (in fiuid dynamics language), the field $ appears 
to be slaved to the field 6. 
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Appendix A. Parameters for DNA modelling 

Our model Lagrangian (1.7) depends on several parameters. Whatever the interest of the 
model per se, and albeit we conducted our discussion in general terms, we are primarily 
motivated by its application to describe the nonlinear dynamics of DNA; hence we are 
interested in the values of these parameters for the physical situation (also to be used in 
numerical plots later on). 

We will not discuss the derivation of these, referring to [5, 24, 29, 31, 34] as well as to 
the detailed discussion given in [4]. 

As for the geometrical parameters, the physical values are 

A ~ 10. 5A , R ~ 6.0l , r ~ 3.0l . (^.1) 

Note that io = 2{A — R — r) is the length of the pairing H-bonds at equilibrium; we have 
^0 — 3.0 A, in agreement with the physical situation [31]. 

In our discussion of the model, we have chosen to adopt Yakushevich contact approxi- 
mation io = 0, which entails A = R + r. We choose then to rescale ^ by a factor (5/6), 
so to satisfy this. Note that as our formulas only involve r and R, and not A, this does 
not show up in there. 

The dynamic parameters are the effective masses M and m; their values are related 
with the moments of inertia Ig and /;, for the rotations considered in the model. For these 
moments of inertia we adopt the values [4] Ig — 2 ■ 10~^^ gA and I^ ~ 5.5 • 10~^^ gA , 
which yields for the effective masses the values 

M = /^/i?2 = 5.5-10"^°g , m = /fe/r^ ~ 6.0 • 10^2° g . {A.2) 

Finally, for the coupling constants we adopt (see again the discussion in [4]) the values 

Kp ~ 0.022eV/A^ , Kpr'^ ~ 0.2eV , Kg ~ 0.25eV . {A.3) 

These values were used in the numerical computations and in the plots of solution 
solutions given in the main text. 

In particular, with these values of the physical parameters and for v = 0, the parameters 
K = 2R^J Kpl\J\ controlling the soliton's shape, and a appearing in (B.l) below take the 
value 

K ~ 0.53 ; a ~ 0.35 . (A.4) 

Appendix B. Order e corrections for "d 

In the main text we provided the leading order terms for the solution to (3.4), i.e. a 
solution at order e" for t?, and at order e for 99; it is appropriate to mention that the 0(e) 
correction to the i? computed above, call it er/, can be explicitly computed. 
Indeed the 0{e) terms in the first of (3.4) provides, using (3.9), 

T]" = — a cosT?o V ~ P sini?o ; (^-1) 

we have written, for ease of notation, 

a = 4Kp/{fi + Mv^) , l3 = SKpMrov^ /[A{fi + Mv^f] . 
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We were not able to solve the equation (B.l) in general; however, if we look at the case 
V ^ 0, then /3 ^ and the equation reduces to 

rj" = ^cos??o??. {B.2) 

It is convenient to express r/(^) as a function /o('i?o) of t?o (which is a monotone function 
of ^, as seen above); with this, and with some simple algebra, the equation (B.2) reads 

2(l-cos'(?o)p" - (sini?o)p' + icos'do)p = 0. {B.3) 

Note that we should as well require that p{0) = = p{2tt), as the boundary conditions 
are satisfied by the first order contribution iDq. 

The solutions to (B.3) satisfying the boundary conditions are given in terms of the 
hypergeometric function 2^^! by 

pi^o) = V(l-cos^o) 2Fifi^^,i^ti^,l,^(l-cos^o)j ; (BA) 

this is readily converted into ?7(^) using (3.9) and recalling that 

viO ■■= P[M0] ■ iB.5) 

Note it follows from (3.9) and trigonometric manipulations that 

\ — Qx -\- X 

cos(t?o) = — 7 KTK — , where x := exp[i^^l . (B.6) 

(1 + x^)^ 

As for the stability of this solutions, it follows from our selection of the function space 
to which solutions must belong, i.e. from the boundary conditions at ±00. In fact, we note 
that as T?o(?) already satisfies the appropriate limit conditions for ^ ^ ±00, the correction 
ri{S^) should vanish for ^ -^ ibcx); this implies the stability of the solution. 

The behavior of solutions for large |^| can also be considered as follows. For ^ ^ ±00 
the term cos(i?o) goes to 1, and asymptotically in ^ eq.(B.2) reduces to 

// ^Kp 

\P\ 

where we have used p < 0; needless to say, the solution of this is given (in the ^ -^ ±00 
region) by 



f] = A±smhi2^Kp/\p\0 + B±cosH2^ Kp/\p\0 ■ 
The limit conditions require that A± = zizB±. 

Appendix C. Perturbative expansion for general solutions 

We have considered expansion around Yakushevich soliton within the function space J- 
identified by the boundary conditions (2.4) and by the travelling wave ansatz (2.1). It 
is of course possible to consider expansion around the Yakushevich soliton 'dQ{x — vt) 
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also for functions which are in T but are not necessarily travelling waves. In this case 
the expansion (3.1) is replaced by (we could as well insert a term ^(x,t) of order one in 
$(x, t); the vanishing of it would however obviously result from the equations obtained at 
order e'^) 

$(x,t) = e4)(x,t) + 0(e2) , 

G(x,t) = -dQix-vt) + er]{x,t) + 0{e'^) . 



,2, (C.l) 



Note that (2.4) imply that the perturbation fields 0(x, t) and ri{x, t) must vanish for 
X -^ iboo. 

Inserting the (C.l) into the Euler-Lagrange equations (1.13) we obtain at order e two 
equations. One of these is actually an algebraic relation between (j) and t?0; which reads 
explicitly 

AAKpMr^v^ 
Kc\\M + m)f ^ — Kgb'^) 

We stress that this imply that the $ field is slaved, even without the restriction to travelling 
waves. 

As for the other 0{e) equation, this reads 

iK,5')v,. - {M + m)^u = 4Kp cos(^o)r/ + ^ ^(^^m + m)v'^ - K,5^) '^^^^"^ " ^^'^^ 

Here and above we have used the writing sin(T?o)) cos(t?o) albeit these are known functions 
of ^ = (x — ft), see (B.6), for ease of writing. 

For large \t\ and finite x the (C.3) reduces to the autonomous equation 

{Ks6^)Vxx - {M + m)r]u = AKp cos(i?o)?? ; {CA) 

the dispersion relations for this are easily obtained setting ri{x,t) = fku)exp[i{qx — ujt)] 
and are given by 

V M+m ^ ' 
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